#include "random.h"

#include <ctime>

Random::Random()
{
	init(time(0));
}

Random::Random(uint seed)
{
	init(seed);
}

uint Random::randUInt()
{
	if(index >= n)
		update_state();

	uint result = state[index++];
	
	// Temper
	result ^= (result >> 11);
	result ^= (result << 7) & 0x9d2c5680UL;
	result ^= (result << 15) & 0xefc60000UL;
	result ^= (result >> 18);

	return result;
}

uint Random::randUInt(uint low, uint high)
{
	assert(low < high);

	return low + randUInt() % (high - low);
}

int Random::randInt()
{
	return int(randUInt());
}

int Random::randInt(int low, int high)
{
	assert(low < high);

	return low + randUInt() % (high - low);
}

float Random::randFloat(float low, float high)
{
	assert(low < high);

	// Random float in range [0, 1)
	float r = float(randUInt()) / (float(max_uint) + 1.0f);
	
	return low + r * (high - low);
}

// Box-muller transform based normal distribution
float Random::normal()
{
	static float cached = INFINITY;
	float result, r, t;

	if(cached != INFINITY)
	{
		result = cached;
		cached = INFINITY;
		return result;
	}

	r = sqrt(-2.0f * log(randFloat()));
	t = 2.0f * pi * randFloat();
	result = r * cos(t);
	cached = r * sin(t);
	return result;
}

float Random::normal(float mean, float variance)
{
	assert(variance > 0.0f);

	return sqrt(variance) * normal() + mean;
}

void Random::init(uint seed)
{
	state[0] = seed;
	for(index = 1; index < n; ++index)
	{
		uint t = (state[index-1] ^ (state[index-1] >> 30));
		state[index] = 1812433253UL * (t + index);
	}
}

// TODO: investigate if it's possible to do the same with 1 loop
void Random::update_state()
{
	int kk, y;
	
	const uint upper = 1UL << 31;
	const uint lower = ~upper;
	const uint matrix = 0x9908b0dfUL;

	for(kk = 0; kk < n - m; ++kk)
	{
		y = (state[kk] & upper) | (state[kk+1] & lower);
		state[kk] = state[kk + m] ^ (y >> 1); 
		state[kk] ^= (y & 1) ? 0 : matrix;
	}

	for(; kk < n - 1; ++kk)
	{
		y = (state[kk] & upper) | (state[kk+1] & lower);
		state[kk] = state[kk + (m - n)] ^ (y >> 1);
		state[kk] ^= (y & 1) ? 0 : matrix;
	}
	
	y = (state[kk] & upper) | (state[0] & lower);
	state[kk] = state[kk + (m - n)] ^ (y >> 1);
	state[kk] ^= (y & 1) ? 0 : matrix;

	index = 0;
}

